Thermodynamic properties of a small superconducting grain 
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The reduced BCS Hamiltonian for a metallic grain with a finite number of electrons is considered. 
The crossover between the ultrasmall regime, in which the level spacing, d, is larger than the 
bulk superconducting gap, A, and the small regime, where A > d, is investigated analytically 
and numerically. The condensation energy, spin magnetization and tunneling peak spectrum are 
calculated analytically in the ultrasmall regime, using an approximation controlled by 1/lnjV as 
small parameter, where N is the number of interacting electron pairs. The condensation energy in 
this regime is perturbative in the coupling constant A, and is proportional to dN\ 2 — \ 2 lud- We 
find that also in a large regime with A > d, in which pairing correlations are already rather well 
developed, the perturbative part of the condensation energy is larger than the singular, BCS, part. 
The condition for the condensation energy to be well approximated by the BCS result is found 
to be roughly A > ydkJE- We show how the condensation energy can, in principle, be extracted 
from a measurement of the spin magnetization curve, and find a re-entrant susceptibility at zero 
temperature as a function of magnetic field, which can serve as a sensitive probe for the existence 
of superconducting correlations in ultrasmall grains. Numerical results are presented which suggest 
that in the large N limit the 1/N correction to the BCS result for the condensation energy is larger 
than A. 



I. INTRODUCTION AND SUMMARY OF 
RESULTS 

In the macroscopic limit, a system described by the 
reduced BCS Hamiltonian is well treated by the mean 
field BCS method §. When the size of a superconduct- 
ing sample becomes small, two related questions can be 
asked: what is the lower size limit for which supercon- 
ducting properties are observable, and what is the lower 
size limit for the validity of the BCS theory? 

In 1959 Anderson considered the first question, and 
argued that "superconductivity would no longer be pos- 
sible" once the electron spectrum's mean level spacing, 
d, becomes larger than the bulk superconducting gap, A. 
(l/d = A/"(0), the density of states per spin species near 
the Fermi energy, hence d oc 1/Volume.) This statement 
sets a lower limit for the size above which a grain still ex- 
hibits superconducting properties, but at the same time 
states that such a grain can well be much smaller than 
the superconducting coherence length. Superconductors 
in the regime where the level spacing is comparable to the 
gap energy have been studied for many years both theo- 
retically (e.g., Ref. S) and experimentally (e.g., Ref. |Q, 
see also review by Perenboom et. al. ||). 

Recently, Ralph, Black and Tinkham (RBT) per- 
formed measurements on single superconducting nm- 
scale grains in the regimes of A > d and A < d || . These 
experiments, and the considerable amount of theoretical 
work they initiated |7|-[l7||, found various properties in- 
dicative of strong superconducting pairing correlations in 
grains with A > d (to be called "small grains"), but not 
in grains with A < d (to be called "ultrasmall grains"), 
thus supporting Anderson's criterion. These properties 
include (i) a parity dependent gap in the excitation spec- 
trum (the gap exists only for grains with an even number 



of electrons), which is driven to zero by magnetic field 
HU|!;|l4|] ; (h) a difference of order A in the ground state 
energies of even and odd grains jnj[l3) ; and (iii) a first or- 
der paramagnetic transition induced by a magnetic field 
§0- 

Though w^rasmall grains with A < d do not have as 
strongly developed signatures of pairing correlations as 
the small grains with A > d mentioned above, pairing 
correlations nevertheless do exist in such grains, albeit in 
the form of weaker fluctuations, and they can affect var- 
ious physical quantities. For example, Di Lorenzo et. al. 
p8| find that pairing correlations affect the temperature 
dependence of the spin susceptibility of grains also in the 
ultrasmall regime. 

The crossover regime between small and ultrasmall 
grains has also been studied in some detail numerically, 
using a simple reduced BCS model with a discrete set of 
single-particle levels. ]l2| , p^| , p^l7| ] . In particular, it was 
found that the condensation energy, E conc i (i.e. the en- 
ergy gain of the exact ground state relative to the uncor- 
rected Fermi ground state), smoothly crosses over from 
being extensive (proportional to the size of the system) 
for A > d to being intensive for A < d. 

One of the goals of the present paper is to obtain 
further insights into the crossover from the ultrasmall 
regime, which can be treated perturbatively in the di- 
mensionless coupling (A) of the said reduced BCS model, 
to the small regime, which can not. Our point of depar- 
ture is an exact solution, due to Richardson [Is|j20(| , of 
the reduced BCS model of present interest. By analyzing 
Richardson's solution both analytically and numerically 
in the crossover regime, we elucidate in detail when and 
how perturbation theory in A breaks down, how the an- 
swer depends on the system size, and how the standard 
BCS results are recovered in the bulk limit d <§C A. 
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The bulk regime is of course well-known to require a 
nonperturbative treatment; indeed, the BCS result for 
the condensation energy, 

Eg c J = A 2 /(2d) , (1) 

is not analytical in A as A — * 0, since the bulk gap is 
given by 

A(A) = cu D / sinh(l/A) [~2uj D e- 1,x for A « 1] , 

(2) 

where lud = Nd is the bandwidth about the Fermi en- 
ergy within which the pairing interaction acts (typically 
the Debye frequency). This nonanalyticity arises because 
BCS consider the thermodynamic limit of infinite system 
size (N — > oo, d — > 0, at fixed ujj). 

We shall argue that if instead one considers a system 
with a finite number of pairs, say N, the condensation 
energy E con d(X) is an analytical function about A = 0, 
with a finite radius of convergence given approximately 
by A* = 1/lniV. For A < A*(l — A*), corresponding to 
A < d (by Eq. (§) Q). E cond (X) is found to be well 
approximated by the perturbative result, 

^T d (A)=ln2.AW (3) 

On the other hand, the BCS mean-field result, -E^nd °f 
(0), is found to become reliable only for A > 2A*, cor- 
responding to roughly A > \fuTod. Thus, we identify a 
substantial intermediate regime, 

A* < A < 2A* , i.e. d < A < ^n~d , (4) 

in which neither the perturbative result nor the BCS 
mean-field result adequately reproduces E con d (though, 
roughly speaking, the sum E p c ^ d + E^cj does). 

The existence of this intermediate regime, which to 
the best of our knowledge has not been identified be- 
fore, implies that the regime of validity of the BCS 
mean-field approach for calculating E con d is significantly 
smaller than realized hitherto: the crossover level spacing 
(d > A 2 /ujd) beyond which it becomes inadequate is con- 
siderably smaller than the scale (d > A) beyond which 
the BCS approach formally breaks down (in the sense 
of yielding no non-trivial solution to the self-consistency 
equation ||), and up to which strong signatures for 
pairing correlations can still be observed, as mentioned 
above. 

We are also able to pinpoint the reason for the failure 
of the BCS approach in the intermediate regime (||): we 
shall show in detail that E^ d incorporates only contri- 
butions to E con d from the strongly pair-correlated, "con- 
densed" levels within A of the Fermi energy Ep, but ne- 
glects contributions from all the remaining, "weakly pair- 
fluctuating" levels, which extend to a distance top, from 
Ep- Although the latter levels are so weakly correlated 
that their contribution can be calculated perturbativcly, 



essentially yielding E^ d , this contribution turns out to 
be larger than £^„f as long as A < Xy/uJod, and is not 
negligible compared to the E^ d in the whole interme- 
diate regime (0) . [Note though, that E v ^ d would largely 
cancel out when one considers energy differences between 
eigenstates that differ only in the specific placement of a 
small number of electrons in levels near Ep. An example 
would be the ground state energy difference between an 
even and odd superconducting grain, for which the BCS 
approach would be adequate in the intermediate regime.] 

It should be mentioned here that the question of how 
to recover the BCS gap equation from Richardson's ex- 
act solution has been solved by Richardson himself [E2| , 
by effectively doing a 1/N expansion around the bulk, 
thermodynamic limit. Our work differs from his in that 
we do an expansion in A around the ultrasmall limit for a 
system of finite size, with A < 1/ In TV as small parameter. 

Using the insights gained from our studies of the con- 
densation energy, we also calculate various other thermo- 
dynamic properties of ultrasmall grains at zero temper- 
ature, using a controlled analytical approximation with 
A < 1/lniV as small parameter. Specifically, we cal- 
culate the spin magnetization and susceptibility curves 
and tunneling peak spectrum of ultrasmall grains, and 
find that pairing correlations have their signature in all 
the above physical quantities, even in the regime A < A* 
where pairing correlations are weakest. 

The condensation energy can, in principle, be mea- 
sured by integrating the spin magnetization as a func- 
tion of magnetic field (H) and comparing it to the lin- 
ear curve of a normal grain. In fact, as we discuss in 
Sec. |ll], since the energy levels in the grain are not 
equally (or systematically) spaced, one needs to do the 
measurement on an ensemble of grains. Calculating the 
spin susceptibility of an ultrasmall grain, we find that for 
H 3> d/fiB, pairing fluctuations of levels far away from 
Ep result in a correction of order \ 2 d/ ' [ibH to the nor- 
mal susceptibility. Interestingly, this correction persists 
for all fields H < uid/hb, i- e - well beyond the Clogston- 
Chandrashekar field HbHcc = A/V2 @, at which, for 
bulk systems, a first order transition occurs from the su- 
perconducting ground state to a paramagnetic ground 
state. (Only for H > lod/^b, the grain becomes effec- 
tively "normal" , since then all the levels within u>p> from 
the Fermi energy become unpaired.) The correction to 
the spin susceptibility results in a re-entrant behavior of 
the differential susceptibility as a function of magnetic 
field, which could possibly serve as a sensitive probe to 
detect superconducting correlations in ultrasmall grains 
p8[ . (The consequences of pairing correlations in the 
regime H > Hqc have also been studied by Aleiner and 
Altshuler pi] , who found an anomaly in the tunneling 
density of states). 

Similarly, we argue below that in ultrasmall supercon- 
ducting grains, pairing fluctuations involving levels far 
away from Ep are sufficiently strong that they also leave 
their mark in the specific heat (even for T3>T C ), and in 
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the tunneling peak spectrum. 

All our calculations are done for grains with an even 
number of electrons. The results for grains with an odd 
number of electrons are similar in the ultrasmall regime, 
and will be discussed shortly for each calculated quantity. 

The paper is arranged as follows: 

In Sec. H we calculate the condensation energy of an 
ultrasmall superconducting grain in the regime A < d, 
and also analyze the intermediate regime of Eq. (|) for 
larger grains. In Sec. Ill the spin magnetization of ultra- 
small grains as function of magnetic field is calculated. 
It is shown that the condensation energy is given by in- 
tegrating the magnetization from H = to uop,/ j.ip- In 
Sec. we calculate the differential spin susceptibility 
of ultrasmall grains as a function of magnetic field, and 
find that it exhibits a re-entrant behavior. In Sec. |v| the 
tunneling peak spectrum of an ultrasmall superconduct- 
ing grain is calculated. In Sec. VI we present numerical 
results for the contribution of the "condensed" and "fluc- 
tuating" levels to the condensation energy. 

The technical aspects of our calculations are collected 
in three appendices. In App. [A| a detailed derivation of 
the accuracy of the condensation energy approximation is 
given. In App. |B]the functional behavior of the prefactors 
of the series expansion of the approximate condensation 
energy is analyzed. In App. ^ the series expansion of the 
exact condensation energy is discussed. 



II. CONDENSATION ENERGY OF AN 
ULTRASMALL GRAIN 

A. Richardson's Equations 



solved exactly. They define for each single particle eigen- 



state of Hq the operator 



This 



operator, for any j, is a constant of motion of the Hamil- 
tonian and takes the value ±1 if the level is singly 
occupied, and otherwise. The many body eigenstates 
of (|^) can therefore be classified into different subspaces 
according to their value of the ^'s, i.e., according to the 
configuration of levels within / which are occupied by 
one electron only. The many body eigenstates and the 
eigenenergies of (||) are then found separately for 
each of the above subspaces. 

The electrons in the singly occupied levels are not scat- 
tered to other levels by the interaction term, and the 
singly occupied levels are "blocked" to pair scattering, 
and we therefore designate them by B. The dynamics of 
the singly occupied levels is also trivially given by Hq. 
Therefore, for each set B one has to solve the reduced 
Hamiltonian: 



H 



u u 

3 i,j 



(6) 



Here 6j = cj + cj_ creates a pair of electrons in level j, and 
U is the set of paired levels within /, i.e. the set of all 
levels that belong to / but not to B (the notation, in gen- 
eral, follows Ref. Below, sums over levels labeled 
by j are to be understood as sums over levels within U. 

Once the configuration of unpaired electrons is given, 
Richardson and Sherman ( ^0|, see JT7J] for a review) 
show that the eigenstates of the system are given by: 

i«> = II c Ui**>. \*k)=cnt=i B t\o), 

ieB 



We consider the reduced BCS Hamiltonian: 

E 'J ( W- Xll ^2< : ,.<-, '■, (5) 

J,o"=± i,j 

for a grain with a given, finite number of electrons N. 
The first term is the kinetic term, which we will refer to as 
Hq, and the second term is the interaction Hamiltonian, 
denoted Hj. The sum in Hj is over all the levels inside 
the range Ep — ujp, < e < Ep + u>p>, which we designate 
by /. The Hamiltonian (j|) is the usual BCS Hamilto- 
nian used when discussing superconducting grai ns M |TJ] 
and its validity is discussed in, e.g., Refs. |l^j2j,|2^T (In 
particular, for the model to be valid the grain's dimen- 
sionless conductance, g, must be much larger than one). 
In all cases discussed below we consider states in which 
all levels below Ep — uipi are doubly occupied, while all 
levels above Ep + ujp, are empty. Since the dynamics of 
electrons occupying levels outside the range / and their 
contribution to the total energy are trivially given by Hq, 
we will not consider them henceforth. 

Richardson and Sherman |l9|,^0| showed that this 
Hamiltonian, with a finite number of electrons, can be 



4 = E 



2ej — E v 



(7) 



where 2k is the number of electrons occupying the un- 
blocked levels, and |0) is the state with all the levels 
below Ep — ujp, fully occupied and all the levels above 
Ep — LUp, empty (in our model |0) is the vacuum state). 
The energy parameters E v (with v = 1, . . . , k) are the 
solutions of a set of k coupled nonlinear equations, the 
i^'th equation of which is given by: 



1 

Xd 



E 



E,, — E v 



E 



i 



2e,j — E v 



0. 



(8) 



The total energy of the system is given by |19|;|20|] : 

B k 



S = E^+E^ 



(9) 



i/=i 



Since the ground state of a grain with an even number of 
electrons does not contain any singly-occupied levels (i.e. 
U = I), the even ground state energy is simply E g s , = 
£Li E v . Its A -> limit is E g . s .{\ = 0) = £* =1 2e u , 
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where {2e u , v = 1, . . . , k} is the set of the k lowest-lying 
single-pair energies. [This is consistent with the observa- 
tion, following from Eq. (@), that in the limit A — > the 
set of EvS reduces to a set of k single-pair energies 2ej, 
which, for the ground state, must have the lowest total 
energy possible.] Consequently, the interaction energy of 
the even ground state, Ei n t(X), defined to be the reduc- 
tion of the exact ground state energy as the interaction 
is turned on from zero to some finite A, can be written 



E mt (X) = E g . s .{0) - E g . s .(X) 



E 



8E U 



(10) 



where we introduced the energy differences 5E V = 2e v — 
E u . A closely related quantity is the condensation energy, 
E CO nd{X), defined to be the energy gain of the exact even 
ground state relative to the uncorrelated Fermi ground 
state: 

E C ond{X) = EF.g.s(X) — E g , s (X) (11) 

fe 

= j^(2e„ - A - = E int (X) - kXd . (12) 
i/=i 

The A-contribution in the first sum in Eq. (|l^) is the 
Hartree self-energy of level j in the Fermi ground state. 



B. Perturbative results for E con d and Ei, lt 

Let us now consider the case that the set I of inter- 
acting levels consists of 2N equally spaced energy levels 
between Ep — ujd and Ep + u>d, occupied by 2N elec- 
trons, so that k = N . Measuring the single-particle en- 
ergies w.r.t. to the bottom of the interacting band, we 
thus take Cj = jd, where j = 1, . . . 2N and d — ljd/N. 
(Note that N ^ N the total number of electrons in the 
grain, which is of order 2Ep/d, not 2u)p>/d). 

Using Eq. (|J), the energy differences 8E V occurring in 
Eq. (no) can be rewritten as 



where 



SE V = 2e„ 



2N 

E 



Xd 



1 



1 — A a„ 



N 



2e, — E v 



E 



En — E V 



For small A, it is natural to approximate 5E U by 

A 



6E» 



X,,d 



where A„ 



1-AoS 



<r; — a v {X — 0) is given by 

2N , N 



(I 



E 



1 

2j - 2v 



E 

/x=l(#y) 



2[i-2v' 



(13) 



(14) 



(15) 



(16) 



The accuracy of this approximation is studied in App. ^ 
(by deriving an expression for 8a ll — a u — aP v ), where we 
find that the relative error in 8E V depends on both A and 
N: Specifically, we find that for all i/, 

8E v /8El = l + O(l/(lnA0 2 ) for A<l/(21nA), (17a) 
5E V /5E° = l + C(l/c 2 ) for 

1/(2 In JV) < A < 1/lniV- c/(lniV) 2 , (17b) 



for any c > 1. Note that Eq. (|17bJ) implies the emergence 
of a second scale near A = 1/lnA, namely l/(ln./V) 2 . 

To the accuracy given by Eqs. (|l7|), the interaction and 
condensation energies can be approximated by 



E int = E cond + N ^ d 



N 



N 



^5E° = ^\ v d, 



(18) 



where X v is given in Eq. (|T^). This result coincides with 
that obtained by Matveev and Larkin (Eq. (17) of |flo|j); 
moreover, our approach allows us to give a controlled 
estimate of the error introduced by this approximation, 
both for Eq. ([18]) and for our explicit calculation of E int 
and E con d in App. [b|. Interestingly, Eq. (|l^) can be in- 
terpreted as a sum over the Hartree selfenergies X v d of 
the lowest TV levels, each of which is evaluated using its 
own, level-specific "renormalized coupling constant" A„ 
(thus motivating our choice of notation) . The emergence 
of such renormalized coupling constants has been noted 
before j2(| , in particular by Matveev and Larkin [juj and 
Berger and Halperin ]13]j , Matveev and Larkin, for ex- 
ample, were concerned with perturbatively calculating 
a certain parity parameter that was essentially equal to 
Xfjd/2, and found 



A 



ML 
N 



X 



l-AlnCwu/d)' 



(19) 



in agreement with our result [Eq. (fL5[)l for An (see 
Eq. ( p36| ) and the statement following it). 

Now, calculating the interaction or condensation ener- 
gies is considerably more involved than calculating the 
parity parameter of Matveev and Larkin, since, in con- 
trast to their calculation, not only one but all N renor- 
malized couplings A„ enter in Eq. ( |i~8| ) for E® nt or E® ond . 
This is a major complication, since their ^-dependence 
turns out to be sufficiently important that it is not pos- 
sible to replace all A„ by a single "effective coupling con- 
stant" . 

Nevertheless, progress can be made by expanding Ef nt 
or E® ond in powers of A and analyzing the convergence 
properties of the resulting series. This is done in App. [b| 
[for E® nt , but here we shall give the results for E® ond , 
which is slightly more convenient, since it lacks the 
Hartree term]. It is found that the convergence radius 
of the power series for E® ond (X) is 



A* = 1/lniV 



(20) 
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The regime of analyticity, 

A<A*, i.e. A < d , 



(21) 



[by Eq. will be called "regime I" below. Within 

regime I, we obtain an analytical expression for E® ond (X) 
as a series in A. We find (see App. |b|) that the se- 
ries for E® ond does not have one parameter which de- 
scribes the ratio between consecutive terms in the se- 
ries. Denoting the m'th term in the power series as 



E 



0(m) 



A while the high powers fulfill the 
A • InN. This results in hav- 



cond ' we show that the low powers fulfill the relation: 

F 0(m+1) , E 0{m) ^ 
cond I cond 

relation: h cond / E cond 
ing two separate scales in A. While the high powers dic- 
tate the convergence radius of the series to be A*, their 
contribution is large only for A > A*(l — A*) (see App.^), 
introducing the aforementioned second scale of l/(ln N) 2 
near A = 1/lniV. As a result, for A < A*(l — A*) (i.e. 
in most of regime I) , E® ond is well approximated by the 
contribution of the low powers, which turn out to corre- 
spond simply to the second order perturbative result (up 
to a relative correction of l/\nN, see App. H): 



E° CO nd * 2ffi(A) = ^ 2 • X 2 u; D [l + 0(1/ In N)} 

for A<A*(1-A*). 



(22) 



This is illustrated in Fig. |l|. Intuitively speaking, this 
contribution can be attributed to pairing fluctuations in- 
volving allthe levels in the range Ep—ujjj < e < Ep+uj£,. 



C. Analysis of the intermediate regime A* < A < 2A* 



Although we are not able to extend the analytical cal- 
culation to the regime of A > A* (i.e. A > d), we are 
able to draw some conclusions about the value of the 
condensation energy in the latter regime: 



First, we note that the perturbative result (22) for the 
condensation energy at A = A* is larger than the BCS 
mean field result (||) as function of A, i.e. E^ nd (X*) » 
E^£ d (X), as long as A < X^/u)p>d. In this regime 
E^nd 'W i s thus also much smaller than the actual con- 
densation energy E con d(X) (since, assuming monotonic- 
ity of E con d(X) as function of A, we have E con d(X) > 
E cZd( X *) for A > A*). In terms of A and N the con- 
dition is In [In 2 • (A*) 2 /2] + 2/A > IniV, which, for large 
N, is roughly A < 2A*. [Note that the exponential de- 
pendence of A on A causes a relatively small change in 
the condition for A (< A* versus < 2A*) to translate into 
a parametric change in the condition for A (< d versus 
X\fijj£,d)\. The E v ^ d contribution in Eq. ( p4| ) becomes 
significantly smaller (by a factor A 2 ) than the E^ ld con- 
tribution only for A > \ZuTod. Thus, we identify an 
intermediate regime, 
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FIG. 1. The condensation energy of a grain with N = 1024, 
in units of level spacing, is plotted as a function of A. The 
solid line is the numerical solution of the exact Richardson 
equations. The dashed line is the second order approxima- 
tion. The dotted line is the BCS approximation. The BCS 
approximation is good for A ^> A* = 1/lniV. In the inset 
the same graph is given for a small range of A and a much 
smaller range for E con d- The value at which the perturbative 
term equals the BCS term tends asymptotically to 2/lniV 
(see text), but here it is somewhat smaller since N is not very 
large. 



d < A < \/uj D d 



A* < A < 2A* 



(23) 



[by Eq. (g)], to be called "regime II", in which the BCS 
mean-field approach is severely inadequate for calculat- 
ing E con d, but which, according to the three properties 
mentioned in the introduction, nevertheless already fea- 
tures strongly-developed pairing correlations. In other 
words, the condition for the adequacy of the BCS mean- 
field approximation (A > 2A*, "regime III") is more re- 
strictive than the condition for the existence of strongly- 
developed pairing correlations (A > A*). Importantly, 
this also means that the BCS mean-field approach be- 
comes inadequate already for much smaller level spac- 
ings, d ps A 2 /ojp,, than those at which it formally breaks 
down (in the sense of yielding no non-trivial solution to 
the self-consistency equation), which occurs for d > A. 

The inadequacy of the BCS approximation in regime 
II stems from the abundance of "fluctuating" levels com- 
pared to "condensed" levels. Each "condensed" level 
within a range A from the Fermi energy contributes ap- 
proximately A/2 to the condensation energy, and having 
A/d such levels gives the BCS term A 2 /2d. Though each 
"fluctuating" level outside this range contributes only 
an amount of order (dX) 2 /d to the condensation energy, 
there are con/d such levels, and for A < X\fu>r)d the total 
contribution X 2 u>d of all fluctuating levels is larger than 
A 2 /2d. This sets an energy scale, \JuJod, which A has to 
exceed before the BCS approximation becomes reliable. 
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The above interpretation of the relative contributions 
of "condensed" and "fluctuating" levels to the total con- 
densation energy is confirmed by a detailed numerical 



VI. 



analysis, see Sec. 

Second, by numerically analyzing Richardson's equa- 
tions ( |I^ , |2tj| , ^7| , p8| , see also a review in ^5|), we find 
that in the regime A > A*, the condensation energy can 
be written as 



Econd (X) = E™${\) + A + a(\)E?Z l d (\) 



pert 



(24) 



where a(A) is a function of A of order unity. (A rather 
similar, but not identical, form was obtained in Eq.(44) 
of from a fit to numerical results for E con d{\) ob- 
tained with the density matrix renormalization group.) 
As will be discussed in more detail in Sec. N% the first 
two terms in Eq. ( pi|) represent the contributions of those 
levels lying within A from Ep (to be called "condensed 
levels"), while the last term is due to the remaining lev- 
els within wo from Ep (to be called "fluctuating levels" ) . 
According to Eq. (pi|), the size- independent correction to 
the BCS result (i.e. the leading order \j~N correction rel- 
ative to the extensive, bulk result) is at least A. 

The numerical analysis carried out in Sec. VI and 
App. also give evidence that Ei nt (and also E conc i) 
is an analytical function on the positive real axis of A 
with a radius of convergence around A = of approxi- 
mately I/lniV. This is in agreement with our analytical 
treatment of the perturbation series in App. [b|,|^. 

The results for the condensation energy of grains with 
an odd number of electrons are similar. In the ground 
state of an odd grain the state at the Fermi level is oc- 
cupied by a single electron. Due to the considerations 
above, one does all the calculation neglecting this level, 
and therefore, when the ground state energy is concerned, 
a grain with an odd number of electrons is equivalent to 
a grain with an even number of electrons with a noninter- 
acting level spectrum which does not contain the single 
level at the Fermi energy, and otherwise identical. This 
change introduces only small quantitative changes in the 
results above. 

One way in which one can, in principle, measure the in- 
teraction energy of an ultrasmall superconducting grain 
is by a measurement of the specific heat. The interaction 
energy is then given by 



[c s (T)-c n (T)]dT. 



(25) 



c n(s) = dE n ( s j jdT ', where E n ^ s \ is the thermal average of 
the energy of a normal (superconducting) grain. While 
in macroscopic samples one obtains the leading order 
(extensive) term of the interaction energy by performing 
the above integral from zero to T c , in ultrasmall grains, 
since the fluctuations involve states in the whole range of 
Ep — ujpi < e < Ep + wj), one has to replace the upper 
limit of the integral by T max w ujp in order to have a 
good estimate of Ei nt . At T > T max , one expects that 
the interaction term in the Hamiltonian would play a 



negligible role, and (E s ) and (E n ) would be roughly the 
same. 

Another way to measure the interaction energy is by 
spin magnetization measurements, as we discuss in the 
next section. 



III. SPIN MAGNETIZATION OF AN 
ULTRASMALL GRAIN 

Since the condensation energy of an ultrasmall grain 
has contributions from all the levels within the range of 
wo, in order to measure it one has to probe all the levels 
within this range. One way to do this is to put an ul- 
trasmall, preferably pancake-shaped grain in a magnetic 
field parallel to the flat direction. One can then neglect 
orbital magnetization, and consider only the Pauli para- 
magnetism psfl . 

The interaction energy can then be obtained by: 



Ej, 



(M n - M s )dH, 



(26) 



where M„( s ) is the magnetization of the normal (su- 
perconducting) grain. This is a general thermodynamic 
identity, relying only on the fact that the electrons further 
than ujp) from Ep are noninteracting, so that M n (H) — 
M S (H) for > Lop,. We now derive this relation 

for ultrasmall superconducting grains, and calculate the 
magnetization of such grains for H 3> d/fj,p. 

We introduce the Zeeman term to the Hamiltonian (|5|) 
changing ej — > ej — afipH (taking the g factor to equal 
2). Each eigenstate of the Hamiltonian (|5|) is also an 
eigenstate of the modified Hamiltonian, with an energy 
Eh = Ep = Q — /ipHini — n±), where n^{n\) is the num- 
ber of levels singly occupied by an electron with a spin 
in (opposite to) the direction of the magnetic field. 

We consider, as above, an ultrasmall grain with an even 
number of electrons and neglect orders of A higher than 
two in the calculations of the eigenstate energies below. 
At T = and zero magnetic field the ground state of the 
system has no broken pairs, meaning there are no bare 
levels occupied with a single electron. Of all the states 
with I broken pairs, the one with the lowest energy will 
be denoted ipi , and its energy Ei . One can show that ipi 
has all the I levels closest to Ep from above and all the 
I levels closest to Ep from below singly occupied, while 
all the other levels are not singly occupied. For H ^ 
all the electrons in the singly occupied levels will have 
their spin in the direction of the magnetic field. In this 
case Ei{H) = Ei(0) - 2l^ B H. For T = and finite H 
the ground state of the system is that ipi with the small- 
est Ei(H) of all Vs. While for a large superconducting 
grain an abrupt transition from I = to I = A/(v / 2d) 
occurs at H — A/ (y^fip) ]23], in an ultrasmall grain the 
number of broken pairs in the ground state increases by 
one at a time as H is increased [H. The magnetic field 
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for which the transition of the ground state from ipi-i 
to tpi occurs is denoted Hi. For Hi < H < Hi + x, tpi is 
the ground state of the grain with ground state energy 
Ei(0) — 21[1bH, and therefore the magnetization equals 
2lfis- The magnetization is a step function in H, with 
equal steps of magnitude 2\ib- One needs only to find 
the values of Hi to get the magnetization curve. The 
above picture is true also for a normal grain. (By nor- 
mal and superconducting grains we mean here similar 
grains, with the same single particle noninteracting spec- 
trum, that differ only by the value of A, which is zero 
for the normal grain and finite for the superconducting 
grain. The relation of the above to a realistic situation is 
discussed below). From its definition as the solution of 



E s ' n {H s ' n ) 



E^(H° /n ) , 



(27) 



(for both pair-correlated or normal grains), Hf is given 
by 



2 A1Sj fff / " = Sf / "(0)-^(0) . 
It follows that 

f^2„ B H^=E^JO)-E s / n (0) 

1=1 



s/r 



(28) 



(29) 



Taking l max — ^d/<1 and subtracting the equation for 
normal grains from that for pair-correlated ones, we find 



£ 2» B (H; - Hf) = ES(0) - E^(0) = E,, 



(30) 



1=1 



where we took Ef (0) = Ef (0), since at energies 
beyond l m axd = lod the pairing interaction is no longer 
operative. But, as can be seen from Fig. ^ (drawn for 
equally spaced normal and pair correlated grains), the 
sum on the left-hand side equals the area between the 
solid and dashed lines, and hence also equals the integral 
in Eq. @. 

Finding Hi amounts to solving Eq. (p7j). We first as- 
sume that the noninteracting energy levels in the grain 
are equally spaced. For a normal grain this equation then 
reduces to: (21 — l)d — 2\ibH = 0, where the first term 
is the extra kinetic energy of the I state compared to the 
I — 1 state, and the second term is its gain in Zeeman 
energy. In an ultrasmall superconducting grain one has 
to add the energy contributions due to Hi to the differ- 
ent ground states. To second order in A, one can show, 
by using either Richardson's equations or perturbation 
theory, that the difference in the interaction energies of 
V>z and ipi-i is 



Xd 



/ N+l-l N+l \ 

2 E Vi+EVi \ X 2 d~Xd + hi(2l,N)X 2 d, 

\j=2/-l j=2l J 



where we define \n(i,j) = X)fc=i V^- Therefore, the 
equation for Hi is 



(21 - l)d + Xd + ln(2Z, iV)A 2 d - 2[i B H = 



(32) 



The above equation is true for all I < up/d, while for 
larger I the interaction term vanishes, and one obtains 
the same equation as for the normal grain. 

The first term in the equation reflects the kinetic en- 
ergy cost of breaking the Vth pair, and is similar to the 
normal grain case. The second term reflects the 'direct' 
(Hartree) energy cost of breaking a pair, coming from the 
diagonal part of the interaction term in the Hamiltonian 
(||). This term is not I dependent, and therefore is not 
reflected in the susceptibility, as we shall see in the next 
section. The third term is the result of the two levels, 
the one I below Ep, and the one I above Ep , becom- 
ing blocked to pairing fluctuations. Its magnitude is a 
decreasing function of I, since as the levels are further 
from Ep their contribution to the pairing fluctuations is 
smaller. This dependence on / is reflected in the suscep- 
tibility. 

In Fig. |^ we plot the magnetization curve for a normal 
grain (A = 0) and a superconducting grain with the same 
equally spaced noninteracting spectrum. 

Using Eq. (^o|), Eq. ( |32|) and the expressions above it, 
one immediately confirms the equality (to second order in 
A) of the interaction energy, as was calculated in App. 
and the integral in Eq. (|26 



16 
M/hb 



12 



(31) 



MsH/d 

FIG. 2. Magnetization curve of normal (solid) and super- 
conducting (dashed) grains with equally spaced noninteract- 
ing spectrum. The width of the rectangles between the curves 
decreases with increasing magnetic field due to the decrease 
of the second order term. The sum of all the areas of the 
rectangles equals E con d/d. 
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So far we considered the idealized case of grains with 
equally spaced energy levels. In order to relate to exper- 
iment, we now relax this assumption. It is not possible 
to see the effects of superconducting correlations on the 
condensation energy by measuring only a single ultra- 
small grain, since the fluctuations of the noninteracting 
energy levels cause larger shifts in the position of the 
Hi's than those induced by the superconducting interac- 
tion. We therefore consider an ensemble of grains with 
the same noninteracting mean level spacing, d, and an 
energy spectrum that obeys GOE statistics. We assume 
that the pairing interaction constant in all the grains is 
the same, given by Ad and calculate the mean spin mag- 
netization of such an ensemble for H ^> d/fis- For each 
grain, the equation ( |32"1) for Hi now becomes: 

Ci + Xd + ln(2Z, N)X 2 d - 2/x B #i = 0. (33) 

£z is the energy difference between the Z'th level above the 
Fermi energy and the Vth level below it in that grain. 
The Hartree term is not affected by level statistics, and 
we neglect the change incurred by the second order term 
due to the effects of level statistics, since this change 
is small compared to its mean value. We approximate 
the second order term in Eq. ( p3| ) by X 2 dln[u;E, / (2[1bH)] 
(replacing I inside the logarithm by its mean value, and 
replacing In by In). We then obtain, for a given magnetic 
field, for each grain, an equation for I, the number of 
broken pairs. It is given by the maximum k that satisfies 
the equation: 

2 f i B H - Xd- X 2 d\n[cj D /(2fi B H)} > ( k . (34) 

The mean value of the magnetization of a grain at a given 
H is therefore: 

M S {H) = 2n%H/d- Xfi B - AVfl Hud/QubH)] . (35) 

The variation around the mean value is given by the vari- 
ation of the number of levels within the energy range 
given by the left side of Eq. (Q) . Since the level statis- 
tics of the grains is given by GOE statistics, the varia- 
tion in the magnetization of one grain is approximately 
5M S (H) = ^BH2fi B H/d]/Tr 2 (see e.g. @). This varia- 
tion is indeed larger than the shift of the mean magne- 
tization compared to that of a normal grain [Eq. (p5|)1, 
but in an ensemble of n grains the variation reduces as 
while the shift in the mean value does not change. 
One can therefore, in principle, measure the interaction 
energy of ultrasmall superconducting grains by measur- 
ing the magnetization of an ensemble of such grains, and 
calculating the integral in Eq. (p6|). While M s is mea- 
sured, M n is given by the straight line starting from the 
origin with a slope equal to the measured ensemble mag- 
netization at H > loo/ hb- The Hartree term in Eq. ( |35| ) 
shifts the magnetization of a pair-correlated grain rela- 
tive to that of a normal grain by a constant, resulting in a 
parallel line not intersecting the origin. Its contribution 
to the integral is trivially Xlud- The second order term 



changes the slope of the magnetization, and introduces 
a nonlinear correction to the normal Pauli susceptibility, 
which we discuss in the next section. 

The consideration of grains with an odd number of 
electrons would lead to similar results in the regime 
HbH 3> d. The magnetization graph for an odd grain 
would be similar to that in Fig. ||, only shifted by one 
unit down and half a unit to the left, not affecting the 
average quantities discussed. 

IV. RE-ENTRANCE OF THE SUSCEPTIBILITY 

Measuring the interaction energy by a magnetization 
measurement might be a difficult task, since it requires 
very high magnetic fields, of the order of ljd/^b- As 
an alternative, we propose here a susceptibility measure- 
ment which would reveal the presence of superconducting 
correlations in ultrasmall grains, and only requires mag- 
netic fields of order d/fj, B - 

Let x s/n (H 7 T) = dM s/n {H 1 T)/dH denote the spin 
susceptibility as function of magnetic field and tempera- 
ture, for a superconducting or normal grain, respectively. 
Di Lorenzo et. al. calculated x s (0,T), finding that 
even for ultrasmall grains it has a minimum at T w d, im- 
plying a re-entrant behavior as function of decreasing T. 
Since this re-entrance differs from the monotonic increase 
expected for the Pauli susceptibility x"(0,T) of normal 
grains, they suggested that it could be a sensitive probe 
to detect superconducting correlations in such grains. 

In this section we discuss an analogous but complemen- 
tary quantity, namely X s (i?, 0). We find that x s (H,0) 
has a maximum at H « d/fis, and decreases as 1/H 
for H > d/fj, B (see Fig. [|. Thus, x' s {H, 0) shows a re- 
entrant behavior in ultrasmall superconducting grains, 
just as x s (0i^) does. Since this again contrasts with 
the Pauli susceptibility X n {H, 0) of normal grains |pl| , p2l , 
measuring x s (H,0) as a function of H could possibly 
serve as a sensitive probe to detect superconducting cor- 
relations in ultrasmall grains. 

For H 3> dj [Lb we use Eq. ( |35| ) and obtain, to first 
order in d/^sH: 

X s = Xo + X 2 hb/H. (36) 

where xo = The susceptibility is a decreasing 

function of H, and the positive 1/ii-correction to the 
normal grain susceptibility, x n (f° r H 3> djji b , to first 
order in dj '/j,bH, one obtains x™ = xo |3l]j3^]) is smaller 
than the leading, normal, term by A 2 d/2/is-ff = A 2 /2l. 

The intuitive reason for why the correction is positive 
is as follows: For a given magnetic field, the magnetiza- 
tion of a pair-correlated grain is, on the average, smaller 
than that of a normal grain (see Fig. ||), because break- 
ing pairs to increase the magnetization costs pairing en- 
ergy. However, since the pairing energy per extra pair de- 
creases the further the pairs involved lie from Ep, the dif- 
ference between the two magnetization curves decreases 
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with increasing H. Consequently, it requires a smaller 
if-increment to break the next pair for a pair-correlated 
than a normal grain, implying a larger susceptibility for 
the former. 

The result in Eq. (^6|) is already sufficient to estab- 
lish the re-entrant behavior of the susceptibility \ s (H, 0), 
since as H is lowered below d and approaches zero, 
X s (H, 0) decreases and approaches zero too, due to level 
repulsion. Precisely at H — the susceptibility x s (H, 0) 
of an odd grain has an additional 5(H)-\ike peak due to 
the contribution of the single, unpaired electron at Ep] 
in fact, for finite T it is the contribution of this unpaired 
electron that is responsible for the re-entrance of x s (0, T) 
predicted in Ref. jl8[. However, for any non-zero H the 
spin of this electron is fully aligned with the magnetic 
field and hence makes no contribution to X s {H > 0,0). 

We now proceed with a calculation of x s , which gives a 
quantitative estimate of the magnitude of the re-entrance 
effect. We consider an ensemble of odd and even grains. 
For a normal grain, \ n (H) is proportional to the proba- 
bility to have a pair of states (I pair), I above and I below 
Ep (denoted —I) separated by Q = ej — e_; = 2/1^7?, and 
is given by [^ljp2| 



X 



2 Ji R{ 



(37) 



where R(x) is the probability of finding two levels a spac- 
ing x apart regardless of the position of the other levels 
§§ (R(x) = 1 - 0(l/x 2 ) for x > 1). We now con- 
sider a single superconducting grain. The l-th pair of 
this grain would contribute, without pairing interaction, 
to the value of the spin susceptibility at H = 0/(2//b). 
However, due to the extra energy cost of breaking a pair, 
which was discussed in the previous section, the I pair 
contributes to the susceptibility at a higher magnetic 
field. This shift is specific to each grain, as it is a func- 
tion of the energies of the noninteracting levels further 
than i from the Fermi energy (the levels closer than I to 
Ep are singly occupied and therefore do not contribute 
to the interaction energy) . While the energy of the I pair 
is arbitrary, and later is taken to satisfy GOE statistics, 
we now make the approximation that the levels further 
than I from Ep are equally spaced with level spacing 
d, the first ones being d apart from the I levels. Due 
to the smallness of the fluctuations in the GOE ensem- 
ble we believe that the above approximation is not only 
proper for H d/fig, but gives fairly good results also 
for H ~ d/fj-B- Under our approximation, which intro- 
duces a modification of Eq. ( |3~i| ) due to the arbitrariness 
of the energies of the I pair, we find that the I pair will 
contribute to x s (H) a t a O-dependent field H(Q) given 
by 



H 



1 



2^i 



N-l 



Ci+Xd+ 1/2 



1 'v^ 1 



=; Ci + j 



Yd 



XhIxq 
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FIG. 3. Spin susceptibility as function of magnetic field at 
T = 0, for A = 0.28, is shown. As H decreases, X s increases, 
until reaching a maximum of 1.02^0 for H m 1.3d//is, im- 
plying a re-entrant behavior. x n (H) (thin solid line) and the 
high field approximation obtained in Eq. ( |36| ) (dashed line) 
are given for comparison. 

Therefore, for an ensemble of ultrasmall superconducting 
grains as considered 

X s (H) = (2^ B /d)P(2fx B H(C)/d) (39) 

|) with Q replaced by f, 



where -ff(C) is given by Eq 
and 



P{2 t i B H{Q/d) = RiC/dWuBdH/dC)- 1 . 



(40) 



(38) 



The result for this calculation with A = 0.28 is given in 
Fig. | 

For larger, but still small grains, where A > d, the spin 
susceptibility is very different at H < A//j,b- However, 
similar calculations j34j to those leading to Eq. ( |3r^ ) show 
that for A 2 / (dfip) < H < ujd/^b one obtains the same 
result as in Eq. ( p6[ ) . The reason essentially is that in this 
regime enough levels are singly occupied, so the energy 
levels involved in the interaction are sufficiently far for 
the perturbative treatment to be valid. 



V. TUNNELING PEAK SPECTRUM 

Superconducting correlations in ultrasmall grains are 
also reflected in their tunneling excitation spectrum. 

For a normal grain, the tunneling peak spectrum is 
simple, consisting of peaks at the single particle excita- 
tion energies of the grain. When the pairing interaction is 
present, the spectrum is much more complex, containing 
peaks at the energies of all the many body states of the 
grain with one electron added or removed. However, for 
a small coupling constant, most of these peaks are small 
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(proportional to A 2 ), and we do not consider them here. 
Instead, we consider only the "primary" peaks, i.e. those 
that survive as A — > and in this limit correspond to the 
tunneling peaks of the normal grain. A non-zero pairing 
interaction reduces their strength by a factor of 1 — A 2 
and shifts their energy. In this section we are concerned 
only with the energy shift of these primary peaks due to 
the pairing interaction. We show that this shift, being a 
decreasing function of energy, causes, for e ^> d, the mean 
spacing between the primary peaks to be smaller than in 
an equivalent normal grain, and therefore introduces a 
positive correction to their density. 

Consider a grain with equally spaced energy levels hav- 
ing an even number (2k) of electrons and ground state 
\4>2k)i and consider the tunneling at positive energies, 
into any eigenstates \4>2k+i) with 2k + 1 electrons. We 
assume that the Coulomb blockade energy is the same 
for the tunneling to all states with 2k + 1 electrons, and 
henceforth neglect it, since we are interested here only in 
energy differences between tunneling peaks. 

To first order in the tunneling Hamiltonian, tunneling 
will occur whenever there is a finite matrix element be- 
tween any state IV'l/c+i) = c s I ^zfc ) ' wnere s is an index 
labeling single noninteracting levels, and any eigenstate 
l^fc+i) with 2k + 1 electrons. 

If the grain is normal, A = 0, then the only relevant 
eigenstates with 2k+l electrons are those in which all lev- 
els up to Ep are filled with two electrons, and one state, s, 
above Ep is occupied by one electron. We define l^/c+i) 
for either a superconducting or a normal grain as the 
lowest energy many body eigenstate of 2k + 1 electrons 
for which the state s is singly occupied. For a normal 
grain IV'ffc+i) and \<p2k+i) are identical. The spectrum 
of tunneling peaks in a normal grain would therefore 
be identical to the noninteracting single-particle energy 
spectrum of the grain. In a similar ultrasmall super- 
conducting grain (A 7^ 0) , pair fluctuations will affect the 
tunneling spectrum in three ways, (i) The primary peaks 
are shifted, and (ii) are reduced in magnitude due to the 
fact that the overlap of the states |"02fe+i) and \4>2k+i) 1S 
smaller than one. (iii) Many small peaks emerge due to 
the small overlap (of order A) between IV^fc+i) and all the 
other many body eigenstates with 2k + 1 electrons which 
are different from \<f>2k+i)- We w ^ n °t consider effects 
(ii) and (iii) here, and proceed with the calculation of 
the mean spacing between the primary peaks in ultra- 
small superconducting grains, as a function of energy. 

The tunneling of an electron into the Z'th level above 
the Fermi energy costs a total energy of 

£(02fc+i) - E(<fig) =ld+ Xd/2 + X 2 dln(N, l)/2 . (41) 

The first term is the kinetic energy contribution. The sec- 
ond term is the Hartree term, and the factor of 1/2 is due 
to the fact that the number of pairs within lu p> below Ep 
changes by, on average, —1/2 when an electron is added 
to the grain (because we have assumed that the band of 
interacting electrons is spaced symmetrically about Ep, 



and Ep shifts upward by half a unit of d when an electron 
is added to the grain). 

Tunneling an electron into the Z'th level also affects the 
interaction energy, which to second order in A is reduced, 
due to blocking level I, by X 2 dln(N, l)/2 ~ A 2 In (N/l)/2. 
(This second order term can be found, using e.g. pertur- 
bation theory, and calculating the difference in the second 
order interaction energy of a Fermi state with 2k elec- 
trons, and the same state up to a single electron added 
to the Z'th level). Both the kinetic energy and the Hartree 
term leave the distance between nearby tunneling ener- 
gies unchanged (d). However, the second order term be- 
comes smaller with increasing energy, and therefore the 
distance between tunneling energies is smaller than that 
of a similar normal grain. This reduction manifests itself 
in the mean spacing of the primary tunneling peaks in 
an ensemble of ultrasmall superconducting grains. One 
can obtain the mean primary peak spacing of such an 
ense mbl e by a similar procedure to the one we took in 
Sec. |h[ |[y] for the spin magnetization and susceptibil- 
ity. Here we obtain the same result in a simpler way. 
We consider a grain with equally spaced energy levels 
as above. The difference between the tunneling energy 
to state \4> 2 k+i) an( l to state | < / > 2&+ 1 ) 1S approximately 
d + A 2 d[ln(7V/(/ + l)) - \n(N/l)}/2 ~ d[l - X 2 /(2l)]. 
The mean density of primary peaks in an ensemble, for 
I l,£>d, including spin degeneracy is therefore given 
by 

*< f > = 5-T^W (42) 

The functional behavior of the primary peak density 
resembles that of the magnetic susceptibility. In both 
cases the correction to the leading term reflects the 
change in interaction energy as levels further from the 
Fermi energy are blocked. 

Similar considerations for negative energies (tunneling 
electrons out of the grain, going from the ground state of 
2k electrons to states of 2k — 1 electrons) will result in 
a similar shift of the tunneling peaks, and the tunneling 
spectrum being symmetric around Ep. 

We now consider shortly the case of the tunneling pro- 
cess (at positive energies) changing a grain with an odd 
number of electrons to a grain with an even number of 
electrons. In the even to odd case described above, the 
change of the number of singly occupied noninteracting 
levels from zero to one induced an upward shift in the 
primary tunneling peak energies. Most primary tunnel- 
ing peaks in the odd to even case correspond to a change 
in the number of singly occupied noninteracting levels 
from one to two, and therefore an upward shift in the 
tunneling peak energies similar to the even to odd case. 
The only exception is the peak closest to the Fermi level, 
which is a result of tunneling into the singly occupied 
level, resulting in the even grain having no singly occu- 
pied levels. This induces a downward shift of this tun- 
neling peak. The functional behavior of the mean peak 
spacing at e 3> d will be similar to Eq. ( [i"2| ) above. 
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We now relax the assumption of the noninteracting 
energy spectrum of the grain being equally spaced, and 
consider the regime < e <C d. Due to the positive shift 
in energy of the tunneling peaks, as a result of the pairing 
interaction, M(e) is smaller for ultrasmall superconduct- 
ing grains compared to its value for similar normal grains 
(this can also be obtained by conservation of the number 
of primary peaks). This can be seen as a first sign for the 
onset of the gap in the density of states in a macroscopic 
sample. 



VI. NUMERICAL ANALYSIS OF 
RICHARDSON'S EQUATIONS 

In this section we analyze the solutions of Richardson's 
equations (^) for the ground state of a grain with equally 
spaced energy levels. 

The above equations for the energy parameters can be 
studied numerically p8| . The energy parameters are in 
principle complex, and one can show that they are ei- 
ther real, or come in complex conjugate pairs. If the 
ground state of the system is considered, it is found |2^] 
that, for A <C 1/ In N all the energy parameters are real, 
monotonically decreasing functions of the coupling. Tak- 
ing iV to be even, we group all the energy parameters E v 
into pairs labeled by an index a, starting from the largest 
two and counting downward: {_Ejv+2-2a, -E/v+i-2a}j with 
a = 1, . . . , N/2. (The case of odd N can be treated analo- 
gously, except that E\ will remain unpaired then.) Each 
pair of energy parameters are real, until, for some A a 
(Fig. 1 of Ref. |2£|), the a'th pair becomes a pair of 
complex conjugate numbers, which they do in order of 
increasing a (i.e. A Q < A a +i). At the transition point, 
A a , the two energy parameters are equal to the value of 
the lower energy parameter at A = 0: 



E 



N+2-2a 



E 



N+l-2a 



2e 



AT+l-2a- 



(43) 



Each energy parameter is an analytic function for A < A Q 
and has a branch point at A a . 

By solving Richardson's equations (Sh for v = N and 
v = N — 1 with the conditions above ( 43) we find that 



Ai = l/(hx/V + ai), 



where IniV = Ylj=i Vi an d < ai < 1. 
It would seem that the interaction energy, 



Ej 



(2e„ - E v ), 



(44) 



(45) 



would be nonanalytical at this point. However, although 
E^,En^i have a branch point at Ax, their sum is an- 
alytical at this point, due to the cancellation of the 
singularities. The analyticity of the sum, as well as 
the result in Eq. (fl4|), are derived |34j in the following 
way: We first define £ and r\ to satisfy the equations: 



E N = 2e^-i + £ + ir) and En-i = 2t N _i + £ — ir], us- 
ing the fact that Em and £V_i are real (in which case 
•q is imaginary) or complex conjugates of each other (in 
which case rj is real). We insert the above definitions in 
Richardson's equations for v = N and v = N — 1, and 
obtain equations for £ and rj 2 as function of A (similar 
equations are given in |p8| , p5| ). We then expand £ and 
rj 2 in a series in SX — (A — Ax), and solve the equations 
for each order of SX separately. Ax is obtained by the 
solution to the zero'th order, and yields Eq. (|44]). Solv- 
ing for the next orders we find that the coefficients, for 
both £ and rj 2 , of (SX) 1 are, up to a factor of order unity, 
(lnJV) 2 \ This results in £ and rj 2 being in fact expanded 
as a series in SX ■ (In N) 2 . This result both reflects a scale 
of l/(lniV) 2 near A = 1/lniV, which we will return to 
later, and shows that the sum En + E^+i can indeed be 
expanded perturbatively in SX. This result also suggests 
that Ei n t(X) is analytical on the positive real axis. The 
only points at which one can suspect nonanalyticities to 
occur are the different A a 's. The reasonable assumption 
that for a > 1 the behavior of Ei nt near A a is similar 
to that near Ax discussed above, i.e. that the sum of 
the singular energy parameters is analytical, leads to the 
analyticity of Ei nt on the positive real axis. This is also 
supported by contour integration in the complex A plane, 
as we discuss in App ^|. 

However, Ei nt is not an analytical function in the whole 
complex A plane. By numerically analyzing Richardson's 
equations as a function of complex A, we find that the in- 
teraction energy has singular points in the complex plane, 
the closest to the origin occurring at approximately 



X smg (N) = l/lnJV± 1.33i/(lniV) 2 



(46) 



In Fig. |j we plot the real and imaginary parts of the 
closest singular points as function of N, with the corre- 
sponding fits. 

This numerical result suggests that the radius of con- 
vergence of the perturbation series for the interaction 
energy around A = is roughly 1/lniV, in agreement 
with our analytical treatment of the perturbation series 
in App. |b]]c| It also reflects a convergence radius of the 
order of l/(ln7V) 2 for the sum En+En-i around A = Ax, 
in agreement with the l/{h\N) 2 scale mentioned above. 

The second part of this section is devoted to estab- 
lishing a connection between the analytical properties of 
Richardson's energy parameters and the BCS theory. In 
particular, we show that in the regime where d <C A 
the BCS result for the condensation energy (|l]) is closely 
related to the singular contribution of the complex en- 
ergy parameters to the condensation energy, and that the 
points, on the positive real axis of A, at which the energy 
parameters become complex, are related to the values of 
A at which additional states become "condensed" (i.e. 
come to lie within A of Ep). 
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0.2 r 

Im A 

0.15 



energy into two parts, one containing the decrease which 
each E2N+2-2a + E2N+i-2a underwent as A is increased 
from to A Q (the perturbative regime), the other contain- 
ing its further decrease for A > A Q (the singular regime). 

Defining, for each energy parameter E Ul its value at the 
branch point as E b v (by Eq. dp) E\ } = E\ r _ A = 2e 2 /-i) 



we write E C C ° Q ^ = 



Tpsmg 
cond 



_ E b 2l = I 
E co% P , where 



cond 



rppert 
comp 



\-E°, 



(48) 



A 

(Re,Im) 

0.3 



0.1 0.2 0.3 0.4 

Re A 

(a) 




(b) 



200 

N 



FIG. 4. (a) The singularity point of the complex coupling 
parameter which is closest to the origin is plotted for (from 
top right to bottom left) N = 16, 24, 32, 48, 64, 96, 128, 192. 
Not plotted is their mirror image across the x axis, with neg- 
ative imaginary parts, (b) The real (top) and imaginary (bot- 
tom) parts of the complex coupling parameter at the above 
singularity points are plotted as function of N. The dashed 
lines are fits, 1/lniV for the real part, and 1.33/ (In N) 2 for 
the imaginary part. The error in both graphs is 0.002 (not 
plotted). 

We separate the energy parameters to two groups, R\ 
and C\, containing those energy parameters which, for 
a given A, are real or complex, respectively. We define 
the contribution of the complex energy parameters to the 
condensation energy (disregarding the Hartree term): 



jpcorap 
cond 



(2e v - A - E v ) 



(47) 



and find that \E™™£ - (A 2 /2d + A)| < 3d for 8 < N < 
1024, for all A < 0.3. The relative correction decreases 
with N and A, and is 0.1 percent for iV = 1024 and 
A = 0.3. 

For those energy parameters which have already be- 
come complex at a given A, we can separate the con- 
tribution of E2N+2-2a + E 2 N 2a to the condensation 



F sing ;o 
cond 



is the singular contribution to E con d or Ei nt . (Note 
that the Hartree term is subtracted in E^^„ when calcu- 
lating the condensation energy; for the analogous calcu- 
lation of the interaction energy, this subtraction should 
be omitted.) Remarkably, our numerical analysis shows 
that this singular contribution is well approximated by 
the BCS expression for the condensation energy (-E^ond = 
A 2 /2d(l + 0{d/A) 2 ) for 64 < N < 1024 for all A < 0.3) 
as is shown in Fig. |(a). This suggests that the BCS 
approximation is equivalent to considering only the sin- 
gular contribution to the condensation energy. By taking 
first the limit TV — ► oo, one would indeed get the singular 
point of the interaction energy to be at the origin (|4q) , 
and no contribution from the perturbative terms. In this 
limit the Hartree term in the reduced Hamiltonian (|^) 
vanishes, and the condensation energy equals the inter- 
action energy. 

These results suggest that the 1/N correction to the 
BCS result in the large N limit is at least A. However, 
one also has to add the contribution of the energy param- 
eters in group R\, which correspond to the levels between 
Ep — lod and Ep — A. This contribution would give an 
additional correction of order X 2 uj]j. Summing the above 
three contributions to E con d, we obtain Eq. (|24|). 

Since E%% * A 2 /2d +A + O(d), and E™« ~ 
A 2 /2d + 0(d), we find that E^ p ~ A [see Fig. |(b)]. 
This implies, since 26^ — E h v equals 2d for even v and 
for odd v, the approximate equation 



2n c (X) 



A{X)/d 
O^A) 



(49) 



where A(A) is given in Eq. (g), and nc{X) is the number 
of pairs of energy parameters {£^v+2-2a, £^+1-20} that 
have already turned complex for the given A. Remark- 
ably, Eq. ( |49| ) tells us that the associated number of bare 
levels £Ar+2-2a and e]v+i-2a i namely 2nc(A), is just the 
number of bare levels within d of A (up to a factor close 
to unity), i.e. the number of what we have called "con- 
densed levels" . The reason for this nomenclature now be- 
comes apparent, since we have just established that the 
singular, BCS-part of E con d arises precisely from those 
energy parameters that have evolved from these 2nc(A) 
bare levels. 
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EPjrt/(2Nd) 




1/ln^ 



0.35 

N 



0.4 
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FIG. 6. The bottom four curves are plots of X a as a func- 
tion of 1/ In (N/a) computed by changing a for fixed (from 
bottom up) N = 128, 256, 512, 1024. The topmost (thick, 
solid) curve is the solution to Eq. d49|). In the inset we 
plot A a as a function of a for (from top left to bottom right) 
N = 64, 128, 256, 512, 1024. The numerical curves (dashed) 
are each fit with the solution to Eq. (^) minus 0.055/ In N 
(solid lines). The error in A a is 0.0002, and is not drawn. 




0.2 

(b) 



0.4 



FIG. 5. The different contributions of the complex energy 
parameters to the condensation energy are plotted, (a) We 
show that -E^ond * s wel ^ approximated by the BCS result, by 
plotting E2Z 9 d /{2N 2 d) for N = 64 (dashed line) and N = 128 
(solid line) and the function exp(—2/\) (dotted line) for com- 
parison. Plots for larger N are not drawn since they are indis- 
tinguishable from exp(—2/X) in the resolution of this figure, 
(b) E^m P /2Nd for N = 64 (dashed line, large steps) and 
N = 256 (solid line, small steps) is plotted and compared 
to exp(—l/X) (smooth solid line), showing that E^m P — A. 
Steps occur at points A a in which pairs of energy parameters 
become complex. 



A 



approx 



(n c ) - X(n c ) ^ 0.055/ In AT, 



(50) 



implying that the approximate prediction becomes very 
accurate for N — > oo. Eq. ( p0| ) represents a one- 
parameter fit that can be used for any combination of N 
and nc, and the quality of which is illustrated in more 
detail in the inset of Fig. [| 

Expanding X a pprox{nc) to first order in \sx{N/nc), we 
find that 



\{n c ) = 1/ In (N/nc) 



1 



lnnc 
hiN ^ (In TV) 2 



(51) 



This equation shows, once again, that the scale of 
l/(lnA') 2 is present also for A > A*, both in BCS theory 
and in the analytical properties of Richardson's equa- 
tions. 



If we solve Eq.(f49|) for A, the result gives a function, 
say Xapprox(nc), that actually depends on N and nc only 
via the ratio N/nc, and that is plotted as function of 
l/]n(N/nc) in Fig. | (thick bold line). This function 
can be used to approximately predict, for given and 
nc, at which value of the coupling constant the nc-th 
pair of energy parameters will become complex. For com- 
parison, we plot in Fig. || also the actual value at which 
this happens, say A(nc)> for several values of A" (non- 
bold lines, obtained by numerically solving Richardson's 
equations). We find that the difference can be well fit by 



VII. CONCLUSIONS 

Though many superconducting properties are limited 
to grains large enough such that d < A, pair correlations 
exist also in ultrasmall grains, where d > A. We calcu- 
lated the effect of pair correlations on the condensation 
energy, spin magnetization and the tunneling spectrum 
of ultrasmall superconducting grains. We found that the 
contribution of pair fluctuations to the condensation en- 
ergy is much larger than the BCS result even for grains 
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in which A > d, and that the condition for the validity of 
the BCS approximation for calculating the condensation 
energy is A > sfuTod. The interaction energy of ultra- 
small grains can, in principle, be experimentally obtained 
through measuring their spin magnetization, which was 
calculated above. The pair correlations result in a posi- 
tive correction to the differential spin susceptibility which 
is proportional to X 2 d/ '(/i^iJ) for /j,bH 3> d, and a posi- 
tive correction to the mean tunneling peak density which 
is proportional to X 2 d/e for e ^> d. The differential spin 
susceptibility at T = of ultrasmall superconducting 
grains shows a re-entrant behavior as a function of H, 
which could serve a sensitive probe for the existence of 
superconducting correlations in such grains. We argued 
that the interaction energy is an analytic function of the 
coupling parameter, with a convergence radius of approx- 
imately 1/lniV. We showed that the BCS result for the 
condensation energy can be obtained from the singular 
part of Richardson's energy parameters, and that the cor- 
rection to the BCS result in the regime where d <C A is 
at least A. 
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Prober for useful discussions. This work was supported 
by the Israel Academy of Science and by the German- 
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APPENDIX A: ACCURACY OF THE ENERGY 
APPROXIMATION 



In this appendix we show that the relative accuracy 
of the approximations in Eqs. (|l5|), ( jl8| ) is as stated in 
Eq. (|l7|). In this appendix we take d = 1. 

From Eqs. (|l6| ) and (14), Sa v = a v — a® is given by: 



2iV 



5a v = - y 



(2j - 2v + 5E v )(2j -2v) 



N 

E 



2{5E ll - 5E y ) 



(2(i -2v + 5E V - 5E fl )(2fi - 2v) 



(Al) 



For A < 1/lnA — c/(ln7V) 2 , we assume (and later 
check for consistency) that < 8E U < 1/c for all 
v. Separating each sum in Eq. (Al) into a sum over 
the levels above and below v : one can see that \8a v \ < 
(1/c) • 3 • 2£ s [2s(2s- l)](- x ) =*> \Sa v \ < 6/c (actually a 
more careful treatment can reduce the numerical factor 
multiplying 1/c to be of order unity, but this is of no 
importance here). 

Therefore, one can write: 



SE U = 



X 



l-Aog + A6„(A)/£ 



(A2) 



where |^(A)| < 6. Manipulating the above equation, we 
obtain 5E V = 5E° + R u where 



X 2 b v {X) 



c-(l-Aa° + A&„(A)/c)-(l-Aa°) 
The relative accuracy is therefore: 

R V {X) . X\b u {X)\ 



r,(A) 



SE V (X) c-(l-Aog) 



(A3) 



(A4) 



Consider now A < 1/lnA - c/(lnA) 2 , and v = N 



(the highest energy parameter). Since 



1 N 



(InN 



InN - l)/2 ~ hxN we see that r N < 6/c 2 . 

One can obtain from Eq. ( filf ) that a° increases mono- 
tonically with v, and by differentiating Eq. (A4) with 



respect to a", that r v increases monotonically with aP u . 
We conclude that r„ < 6/c 2 for all v, and that there- 
fore 5E V = 5E°[1 + C(l/c 2 )] for all v. In the same way 
one finds that 5E V < 1/c for all v, consistent with our 
assumption above. 

By summing over all v 1 one shows that 



Ei, 



£l(l + 0(l/c 2 )), 



(A5) 



where Ef nt is given in Eq. ©. 

To show that the accuracy of the above expression is 
l/(lniV) 2 in the regime A < 1/(2 lniV) one has to assume 
(and later show consistency) that 8E U < 1/ In A for all 
v, and proceed as above. 



APPENDIX B: APPROXIMATE FORMULA FOR 
THE M'TH ORDER TERM OF THE 
INTERACTION ENERGY 



In this appendix we expand Eq. (18) for the approxi 



mate interaction energy E® nt in powers of A, and analyze 
the convergence properties of the resulting series. We 
begin by showing that the order- A m contribution to Ef nt 
has the form 



E, 



0(m) 



dX r 



N 

0\m-l 



d\ m b m N(m-l)\ 
d X m b' m (\n N)( m -V 



for 
for 



m < InN 
m > In ./V 



(Bl) 



(B2) 



where b m and b' m are constants of order unity, and then 
analyze the consequences of this result. Since we are not 
interested here in numerical factors of order unity, we will 
allow ourselves to make some crude approximations. 

The expression obtained (|lj) for the interaction en- 
ergy can be written as: 



/V 



0(m) 



= d£[A 



i/=i 



(B3) 
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where a° is given by Eq. (|l6|). The first order term is First, since the low powers in the series for Ef nt fulfill 
given trivially by E°$ = XdN. 



the relation E^/E^f 



The second order term of the interaction energy is 
given by: 



JV 



i=l 
JV 



1 



N i 

E - 1 



2JV 

^ 2(j - i) ^ ]-i 



2JV 

E 

j=JV+l 



2(j ~ i) 



(B4) 



This result can also be obtained by standard second 
order perturbation theory ( pa ] and references therein). 

For large N one obtains: E°^~k, In 2 A 2 eOV = In 2 A 2 wd. 

The calculation of the higher order terms is more dif- 
ficult. We now make some approximations which enable 
us to find the m'th order within a factor of order unity. 

First we manipulate the a°'s (M) and obtain: 



u-l j N-v ^ IN —v ^ 



(B5) 



_fc=l fe=l fc=JV-i/+l 

This can be approximated by: 

al = - bx[2(N + l)/v - 1] - ln[(jV + l)/v - 1] . (B6) 

For v <C JV one obtains a° ~ — In (N/v)/2, and for 
v ~ N (meaning jV — ^ <C iV) one obtains a° ~ In [N/v). 
We now make a crude approximation: 



JV 



AT 



(l-C-r-^EQnW.y))™- 1 . (B7) 



j'=i 



This approximation is proper only for the v's which 
are either small or close to N. However, since these j/'s 
contribute the most to the sum, we expect this approx- 
imation to be correct within a numerical factor of o rde r 
unity. Indeed, for m = 2, the last sum in Eq. (B7) 



approximately equals N, and the total result we obtain, 
0.5N, is different than the correct result, In 2 • N only 
by a factor of order unity. Increasing to, the approxi- 
mation becomes better, since the relative contribution of 
the levels far from N decreases. For to > 2 we neglect 
the factor (l/2) m in @. For to > In N only the j = 1 



term in Eq. 
by A m • In AT( 



contributes, and is approximated 

For to <C In N one can show that 



JV 



(to- ly.N -\nN ( - rn - 1 '> <J2( ln ( N /j)) m ~ 1 < (m-l)lN 

3 = 1 

(B8) 



and therefore £.L(hi Wi)) m_1 ^ (m - 1)UV. 



Having thus established Eq. (Bl), let us now examine 
its consequences: 



to • A, whereas the high 

0(m) 



powers fulfill the relation E%™ +1) / E^™> ~ A • InJV, the 
series for i?? nt does not have a single parameter describ- 
ing the ratio between consecutive terms in the series. 

Second, while the high powers dictate the convergence 
radius of the series to be roughly 1/lnTV, their contri- 
bution is larger than that of the low powers only for 
A > 1/lnJV- l/(lniV) 2 , introducing a scale of l/(ln7V) 2 
near A = 1/lniV (see also Sec. VI). This can be seen 
by estimating the partial sum of Eq. (B3) for to > In N 
using the result in Eq. B2. Taking b' m = 1 for all to we 
get: 



X m (\nN) r 



(A In TV) 



In JV 



m=ln JV 



lnjV(l- AlnjV) 



(B9) 



For A = 1/lniV — c/(lnjV) 2 , for any c > 1, the above 
sum equals: 



(1 - c/ lnjV) 1: 



(BIO) 



which is smaller than one, while the low orders are pro- 
portional to N. 

Third, for A < 1/ In N- 1/(1iijV) 2 , where the low pow- 
ers dominate, one readily obtains 



E° nt = XNd + E*$ [1 + 0(1/ In JV)] 



(Bll) 



Here the first order term, XNd, is the Hartree term, and 



the second order term, E, 



(2) 



In 2 



J int ~ X ujd, is obtained 

from Eq. (|B4|). The order 1/lnjV can be understood as 
follows. The third order term is smaller than the second 
order term by a factor of order A, which is smaller than 
1/lnjV. All the higher orders are smaller than the sec- 
ond order term by a factor of order l/(ln N) 2 or smaller, 
and since there are about InjV such terms, their sum is 
also approximately 1/ In N smaller than the second order 
term. 

The corresponding perturbative result for the conden- 
sation energy, E® nnd = -E^ond' now immediately follows 
by inserting Eq. ( |B11|) i nto Eq. (Jll|) [with k = N]. The 
result is given by Eq. j22j), namely B^' d (A) ~ In 2-X 2 lod ■ 



APPENDIX C: SERIES EXPANSION OF THE 
INTERACTION ENERGY 

In Eq. ( |l8| ) we obtained an expression for the interac- 
tion energy, whose accuracy for different regimes in the 
range A < 1/lnjV was obtained in App. |a[ We then ex- 
panded the result in a series in A. This series converges, 
in the regime A < 1/lnAT, to the approximate result. We 
were not able to find a series which converges to the ex- 
act result for the interaction energy in the above regime. 
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However, we obtained results which suggest that the in- 
teraction energy is analytical on the positive real axis of 
the coupling parameter, and can be expanded in a se- 
ries in A with a finite convergence radius which equals 
l/(hx/V + b), where b is of order unity (Sec. VI). As 



another check of the above statement we solved, numeri- 
cally, Richardson's equations (||) for complex values of A, 
and calculated the integral: 



E int (z)dz, 



(CI) 



where C is a contour circumventing the positive real axis 
(see Fig. |). 

The integral was calculated for various TV's in the range 
4 < N < 64 (a few contours for each N, each extending 
to a different value of ReX, up to ReX = 0.7). For all N 
the integrals were zero within the numerical error, which 
suggests that there are no singularities on the positive 
real axis. 

The exact interaction energy is given by Eq. @. Ex- 
panding the exact expression as Ei nt = dY^ =1 ct m X m , 
we find that 



(a o)(m-D + j2 b s ( a y 



(C2) 



The approximation we make in Eq. (|18| ) is equivalent 
to taking only the highest power (m — 1) in Eq. (C2) for 
each m, v, and neglecting the sum in the brackets. 



Im X 



Re X 



FIG. 7. An integration contour in the complex A plane. 
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